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SHORT  COMMUNICATION 


Modeling  Loop  Reorganization  Free  Energies  of 
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Solvent  Models 
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ABSTRACT  The  treatment  of  hydration  effects 
in  protein  dynamics  simulations  varies  in  model 
complexity  and  spans  the  range  from  the  computa¬ 
tionally  intensive  microscopic  evaluation  to  simple 
dielectric  screening  of  charge-charge  interactions. 
This  paper  compares  different  solvent  models  ap¬ 
plied  to  the  problem  of  estimating  the  free-energy 
difference  between  two  loop  conformations  in  acetyl¬ 
cholinesterase.  Molecular  dynamics  (MD)  simula¬ 
tions  were  used  to  sample  potential  energy  surfaces 
of  the  two  basins  with  solvent  treated  by  means  of 
explicit  and  implicit  methods.  Implicit  solvent  meth¬ 
ods  studied  include  the  generalized  Born  (GB)  model, 
atomic  solvation  potential  (ASP),  and  the  distance- 
dependent  dieletric  constant.  By  using  the  linear 
response  approximation  (LRA),  the  explicit  solvent 
calculations  determined  a  free-energy  difference 
that  is  in  excellent  agreement  with  the  experimen¬ 
tal  estimate,  while  rescoring  the  protein  conforma¬ 
tions  with  GB  or  the  Poisson  equation  showed  incon¬ 
sistent  and  inferior  results.  While  the  approach  of 
rescoring  conformations  from  explicit  water  simula¬ 
tions  with  implicit  solvent  models  is  popular  among 
many  applications,  it  perturbs  the  energy  landscape 
by  changing  the  solvent  contribution  to  microstates 
without  conformational  relaxation,  thus  leading  to 
non-optimal  solvation  free  energies.  Calculations 
applying  MD  with  a  GB  solvent  model  produced 
results  of  comparable  accuracy  as  observed  with 
LRA,  yet  the  electrostatic  free-energy  terms  were 
significantly  different  due  to  optimization  on  a  po¬ 
tential  energy  surface  favored  by  an  implicit  solvent 
reaction  field.  The  simpler  methods  of  ASP  and  the 
distance-dependent  scaling  of  the  dielectric  con¬ 
stant  both  produced  considerable  distortions  in  the 
protein  internal  free-energy  terms  and  are  conse¬ 
quently  unreliable.  Proteins  2004;57:645-650. 
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INTRODUCTION 

The  population  shift  of  enzymes  from  their  native  to 
functional  non-native  macrostates  during  complex  forma¬ 
tion  is  a  measure  of  substrate  specificity.  The  effect  of 
conformational  selection  is  nicely  illustrated  by  the  X-ray 
crystal  structure  of  the  reaction  of  acetylcholinesterase 
(AChE)  from  Torpedo  californica  with  the  inhibitor  diiso- 
propylphosphorofluor/date  (DFP).1  Structural  changes  take 
place  in  the  acyl  loop  pocket  of  residues  287-290,  showing 
the  backbone  displaced  by  roughly  5  A.  Comparison  of 
experimental  reversible  dissociation  constants  measured 
for  a  variety  of  inhibitors  of  differing  molecular  size 1,3-7 
suggests  that  the  free  energy  penalty  for  loop  displace¬ 
ment  is  on  the  order  of  4  kcal/mol. 

Central  to  the  problem  of  modeling  binding  affinities, 
the  calculation  of  structural  reorganization  is  essential  to 
reliable  and  accurate  predictions  of  complex  forma¬ 
tion.8,9,11  The  experimental  observation  of  significant  loop 
movement  in  the  acyl  pocket  offers  an  excellent  bench¬ 
mark  for  testing  different  methodological  approaches  to 
computational  analysis  of  protein  structures.  This  article 
explores  the  application  of  different  solvent  models  of 
varying  resolution  in  the  calculation  of  the  free  energy 
difference  between  the  conformational  loop  captured  in  the 
DFP-reacted  macrostate  and  native.  Molecular  dynamics 
(MD)  simulations  were  carried  out  with  explicit  treatment 
of  solvent  and  the  free  energy  of  hydration  estimated  by 
application  of  linear  response  approximation  (LRA).  The 
hydrophobic  term  was  calculated  from  a  molecular  surface 
area  (SA)  model.  An  alternative  and  currently  popular 
approach  was  also  applied  of  reevaluating  the  MD  trajec¬ 
tory  by  removing  the  explicit  solvent  and  scoring  the 
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conformations  with  either  a  finite-difference  Poisson  (FDP) 
technique  or  a  generalized  Born  (GB)  model.  A  third 
approach  was  investigated  by  performing  a  MD  calcula¬ 
tion  with  solvent  treated  implicitly  with  a  GB  model. 
Finally,  an  atomic  solvation  potential  (ASP)  and  a  distance- 
dependent  dielectric  constant  model  were  applied  as  fur¬ 
ther  implicit  solvent  schemes  in  MD  simulations. 

While  there  is  extensive  literature  on  GB  and  FDP 
applications,  the  accuracy  of  these  models  still  remains 
uncertain.  An  overall  contrast  of  different  solvent-model¬ 
ing  techniques  is  much  needed  for  an  assessment  of  not 
only  of  the  accuracy  in  predicting  thermodynamic  differ¬ 
ences  between  conformational  states  but  also  modeling  the 
correct  physics  underlying  the  differences. 

METHODS 

Structural  reorganization  between  two  protein  loop  con¬ 
formations  (denoted  as  P  and  P*)  located  in  different 
basins  can  be  described  by  the  reaction 

AG 

(1) 


where  AG  is  the  free  energy  difference  of  the  protein  and 
solvent  system.  The  free  energy  G(0  for  a  macrostate  ^  is12 


Gmacro(Q  =  ~kBT  In  D  =  Gmicro(0  -  TScon{  (0,  (2) 

where  D  is  the  partition  function,  Gmicro  is  the  free  energy 
of  a  single  loop  conformation  taken  to  be  the  most  probable 
microstate  in  0  £conf  is  the  conformational  entropy,  and 
kBT  is  the  Boltzmann  constant  and  absolute  temperature. 

To  determine  Gmicro,  we  partition  the  free  energy  for 
each  conformation  into  separate  contributions;  for  ex¬ 
ample,  conformation  P  gives 


pp 

micro 


=  & 


+  G\ 


[P 

cav? 


(3) 


where  Pint  is  the  internal  energy  of  the  protein,  Gsolv  is  the 
solvation  term  and  Gcav  is  the  cavitation  free  energy  given 
by 


pp 

^cav 


(4) 


where  for  atom  i,  is  the  surface  tension  and  At  is  the 
molecular  surface  area.  The  internal  energy  component  is 
determined  from 

pfnt=pfocal+<dw+wrle,  (5) 

where  Elocal  is  the  bonded  and  torsional  terms  defined  by  a 
given  force  field,  Pvdw  is  van  der  Waals  (vdW)  intra¬ 
protein  interactions,  and  Wele  is  the  energetic  cost  of 
creating  the  charge  distribution  in  an  environment  of  unit 
dielectric. 

Solvent  effects  are  modeled  as  the  sum  of  non-electro¬ 
static  and  electrostatic  interaction  contributions 


^"solv  ^s,ele  ^s,vdW5 


(6) 


determined  for  each  microstate.  From  explicit  solvent 
calculations,  the  free  energy  of  solvation  can  be  estimated 
from  the  LRA  (see,  e.g.,  ref.  9) 

Gfolv  =  a«t/U>P  +  <^ele>F')  +  P(^vdW>P  (?) 

where  U  is  the  potential  energy  surface  of  non-bonded 
interactions  between  the  protein  in  conformation  P  and 
solvent,  {...)P  is  the  ensemble  average  over  charged  P 
microstates,  and  is  the  average  over  uncharged  P', 

and  a  and  (3  are  scaling  constants.  The  form  of  eq.  (7)  can 
be  approximated  by  noting  that  the  solvent  in  the  un¬ 
charged  P'  state  does  not  experience  the  charge  distribu¬ 
tion  of  the  solute,9,10  thus  setting  the  term  (Gf  eie)p'  to  zero. 

Alternatively  to  LRA,  loop  conformations  from  explicit 
solvent  simulations  can  be  rescored  via  either  GB  or  FDP 
after  removal  of  water.  A  computationally  more  tractable 
approach  to  explicit  water  calculations  is  the  implicit 
solvent  schemes  of  MD-GB/SA  and  ASP,  or  simpler  meth¬ 
ods  of  linear  distance-dependent  dielectric  constants. 

The  conformational  entropy  in  macrostate  ^  given  by  eq. 
(2)  can  be  estimated  in  the  quasi-harmonic  approximation 
calculated  from  the  covariance  fluctuation  matrix  C:  13-17 

Cfj  =  <(xf  -  <xf»(xf  -  <xf»)e,  (8) 

where  the  values  of  x  are  atomic  coordinates  of  microstates 
of  P  sampled  from  0  Frequencies  vL  of  the  normal  modes 
are  defined  by  the  eigenvalues  \ 

(2™t)2  =  kBT/\h  (9) 

and  the  entropy  £conf  (0)  is  approximated  by13 
TSC onf(£)  =  y  -  knT  In  [1  -  exp(  -  hv,  /  kBT )] 


+  hvt  /  [1  -  exp(  -  hvi  /  kBT )]  (10) 

where  h  is  Planck’s  constant. 

Calculation  of  native  microstates  started  with  the  Carte¬ 
sian  coordinates  taken  from  the  crystallographic  PDB  file 
2 ACE.  Crystallographic  waters  were  ignored  in  all  calcula¬ 
tions.  The  alternative  loop  conformation  (residues  287- 
290)  was  taken  from  the  2DFP  structure.  To  model  a 
common  structure  differing  only  in  the  loop  conformation, 
the  alternative  loop  was  built  into  2ACE  by  substitution  of 
the  native  and  optimized  by  a  Monte  Carlo  simulation.18 
Figure  1  shows  a  superposition  of  native  loop  and  the 
reorganized  loop.  Backbone  geometries  of  the  two  conforma¬ 
tions  differ  by  the  180°  reorientation  of  the  plane  of  the 
peptide  bond  between  Ile-287  and  Phe-288.  The  peptide 
hydrogen  of  Phe-288  is  pointing  inside  the  active-site 
gorge  for  the  native  conformation  and  outside  the  gorge 
upon  displacement.  Because  of  this  reorientation,  the 
native  loop  is  designated  the  ‘in’  conformation  and  the 
alternative  conformation  as  the  ‘out’  conformation. 

To  determine  A Gmacro  for  ‘in’  — >  ‘out’  loop  reorganization, 
the  computational  approach  builds  upon  a  previous  simu¬ 
lation  study,  which  showed  that  no  dynamic  connectivity 
exists  between  the  two  macrostates.18  Moreover,  both 
states  exhibit  Gaussian-like  probability  distribution  func¬ 
tions  along  the  normal  modes.18  Because  of  the  deep 
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Fig.  1 .  Molecular  illustration  of  loop  reorganization  in  acetylcholinest¬ 
erase  resulting  from  inhibitor  reaction.1  The  ‘in’  loop  conformation  (resi¬ 
dues  287-290)  is  shown  in  blue  and  the  ‘out’  conformation  in  red. 


funnel-like  basins,  the  ensemble  of  protein  configurations 
for  the  ‘out’  loop  with  and  without  DFP  is  nearly  identical, 
and  this  allows  the  application  of  a  two-state  simulation 
model  without  modeling  the  bound  inhibitor.  Without  the 
advantage  of  funnels,  either  structural  constraints  are 
required  on  the  protein  configuration  to  remain  in  the  ‘out’ 
state  or  bound  DFP  must  be  modeled  and  a  computational 
strategy  implemented  similar  to  that  described  by  Warshel 
and  coworkers,8’9  among  others.11  Alternatively,  the  tran¬ 
sition  between  the  two  states  can  be  described  in  terms  of  a 
binding  cycle  of  different  ligands  that  bind  the  native 
conformation  versus  the  non-native.19  The  advantage  of 
this  latter  approach  when  developed  in  a  LRA  framework 
is  the  need  to  only  model  the  change  in  ligand-protein 
interaction  as  a  reasonable  approximation  for  the  corre¬ 
sponding  change  in  the  interaction  between  different 
protein  conformations.19 

MD  simulations  were  performed  with  the  Tinker  pro¬ 
gram20  using  the  AMBER94  force  field.21  The  simulation 
volume  was  constructed  from  a  shell  of  residues  a  within 
12  A  distance  of  residues  287-290  in  the  starting  struc¬ 
tures,  giving  a  total  of  127  active  residues.  The  remaining 
residues  were  tethered  to  their  initial  positions.  Explicit 
water  calculations  of  the  protein  active  region  contained  a 
20  A  shell  of  solvent  (1675  waters)  modeled  by  the  TIP3P 
potential.22  A  restraining  potential  was  used  at  the  solvent- 
vacuum  boundary.  The  reaction  field  beyond  the  boundary 
was  added  to  the  simulation  results  as  a  correction  term 
using  FDP  (see  below).  Similar  protein  regions  were 
simulated  with  a  GB  model,23  ASP24  and  a  modeled  linear 
distance-dependent  protein  dielectric  constant  (e  =  r).  For 
comparison  purposes,  a  vacuum  (e  =  1)  calculation  was 
performed. 

Non-bonded  interactions  were  calculated  using  a  spheri¬ 
cal  cut-off  of  20  A.  The  protein  dielectric  constant  ( ep )  was 
set  to  1,  and  for  the  GB  model  80  was  used  to  model  bulk 
water.  Ionization  states  were  set  corresponding  to  a  neu¬ 
tral  pH.  Starting  velocities  were  assigned  from  a  Boltz¬ 
mann  distribution  at  298K.  The  integration  timestep  was 


set  at  1.0  fs,  and  simulations  were  initiated  with  a  50  ps 
equilibration  phase  followed  by  a  500  ps  production  calcu¬ 
lation,  with  a  total  of  1000  conformations  culled  from  each 
simulation. 

Conformational  entropy  was  determined  by  solving  eq. 
(10),  using  the  solutions  to  the  covariance  fluctuation 
matrix.  All  heavy  atoms  of  the  active-site  region  were  used 
in  the  evaluation  of  TSconf. 

The  LRA  model  was  solved  with  a  and  p  set  to  0.50  and 
0.16,  respectively.9  Cavitation  energies  were  determined 
from  the  molecular  surface  using  the  Connolly  algorithm25 
with  the  solvent  probe  radius  set  at  1.4  A.  The  numerical 
value  of  y  was  defined  as  69  cal  mol-1  A-2  for  all  atom 
types.26  The  vdW  solvation  free  energy  term  for  implicit 
solvent  model  calculations  was  scaled  as  Gsvdw  ~  p>Gcav, 
where  |jl  was  determined  from  LRA. 

The  FDP  method  implemented  was  from  the  program 
DelPhi.27  Electrostatic  potentials  for  each  molecule  were 
calculated  using  molecular  surfaces  to  define  regions  of 
low  dielectric  medium  (ep  =  1)  embedded  in  high  dielectric 
solvent  water  (e^  =  80).  DelPhi  calculations  were  carried 
out  on  a  cubic  grid  of  resolution  0.6  A/point,  giving  a  total 
of  1493  grid  points.  For  the  explicit  solvent  reaction-field 
correction  term,  waters  were  removed  and  the  dielectric 
boundary  approximated  by  scaling  the  atomic  vdW  radii  of 
the  protein  to  encapsulate  the  simulation  volume,  and  a 
grid  of  1.0  A/point  was  applied.  Full  Coulombic  boundary 
conditions  were  applied  for  all  FDP  calculations. 

RESULTS  AND  DISCUSSION 

Table  I  summarizes  the  free  energy  components  of 
AGmacro  determined  for  each  of  the  simulation  models.  The 
reported  errors  reflect  the  statistical  uncertainty  in  the 
sampling  distributions  and  were  evaluated  as  the  stan¬ 
dard  deviation  of  the  free  energy  divided  by  the  square  root 
of  the  sample  size;  namely,  Z(aGmacro)/\/A^,  where  Z  = 
1.6452  is  the  number  of  standard  deviations  from  the 
mean  that  is  needed  to  claim  with  95%  confidence  the 
value  of  the  true  mean.  Convergence  is  an  issue  of  all 
sampling  problems,  and  while  the  simulation  time  and 
sample  size  of  the  calculations  are  moderate,  the  observed 
trends  among  the  solvent  models  have  converged. 

The  result  of  the  explicit  solvent  simulation  using  the 
LRA  to  calculate  the  free  energy  difference  yields  a  A Gmacro 
of  ^5  kcal/mol  and  predicts  stabilization  of  the  ‘in’  state, 
which  is  consistent  with  the  experimental  estimate  of 
A Gexpt  ~  4  kcal/mol. 1,3-7  Because  of  the  complexity  of 
isolating  the  effect  of  conformational  change  and  its  ther¬ 
modynamic  underpinning  from  the  inhibitor  reaction, 
AGexpt  is  a  rough  multipart  estimate.  Nevertheless,  the 
proposal  that  the  AGexpt  is  due  primarily  to  structural 
reorganization  is  corroborated  by  site-directed  mutagen¬ 
esis  of  human  AChE  at  residues  Phe-295  and  Phe-297  in 
the  loop  segment  and  the  effect  on  binding  DFP.6,7  These 
two  phenylalanines  of  AChE  are  conserved  among  species 
and  are  equivalent  to  Phe-288  and  Phe-290  in  Torpedo  (see 
Fig.  1).  Substituting  both  residues  with  smaller  aliphatic 
side  chains  favorably  alters  the  free  energy  of  association 
for  DFP  by  ^4  kcal/mol  and  is  similar  to  binding  inhibitors 
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TABLE  I.  Summary  of  Energy  Differences  (kcal/mol)  from  Simulations 
Employing  Different  Solvent  Models 
for  (in  — >  out)  Loop  Reorganization  on  Acetylcholinesterase* 


AEnn 

AWele 

AG, 

s.ele 

AGsvdw 

AGcav 

-TASmnf 

A  /^fin  — *■  out 
tAUmacro 

Explicit  Solvent 

-1.2 

-116.3 

LRA 

130.9 

0.1 

0.1 

-8.3 

5.3  ±  3.6 

-1.2 

-116.3 

GB 

123.7 

0.1 

0.1 

-8.3 

-1.9  ±  3.2 

-1.2 

-116.3 

FDP 

117.6 

0.1 

0.1 

-8.3 

-8.0  ±  3.4 

GB/SA 

-2.9 

-52.1 

GB 

63.6 

0.1 

0.1 

-4.5 

4.3  ±  3.2 

-2.9 

-52.1 

FDP 

56.0 

0.1 

0.1 

-4.5 

-3.3  ±  3.1 

ASP 

-50.9 

42.8 

ASP 

-19.3 

0.1 

0.1 

-8.9 

-36.1  ±  2.4 

-50.9 

42.8 

GB 

-46.7 

0.1 

0.1 

-8.9 

-63.5  ±  2.7 

-50.9 

42.8 

FDP 

-46.2 

0.1 

0.1 

-8.9 

-63.0  ±  2.6 

e  =  r 

34.1 

-51.4 

£  =  r 

0.0 

0.1 

0.1 

-9.5 

-26.6  ±  2.3 

34.1 

-128.2 

GB 

115.2 

0.1 

0.1 

-9.5 

11.8  ±  2.3 

34.1 

-128.2 

FDP 

101.6 

0.1 

0.1 

-9.5 

-1.8  ±2.4 

£  =  1 

-5.2 

-36.0 

GB 

41.8 

0.1 

0.1 

-14.1 

-13.3  ±2.0 

-5.2 

-36.0 

FDP 

36.7 

0.1 

0.1 

-14.1 

-18.4  ±2.0 

^Experimental  estimate  taken  to  be  on  the  order  of  ~4  kcal/mol.1,3  7  The  term  A Enp  is 
defined  as  bonded  and  torsional  terms  plus  internal  vdW. 


where  no  loop  movement  is  observed  crystallographi- 
cally.6,7  The  simulation  results  of  Table  I  support  the  view 
that  the  energetic  cost  is  consequent  of  loop  displacement. 

Free  energy  terms  A Gcav  and  AGS  vdw  are  insignificant 
in  discrimination  between  the  two  states,  and  A Eexp 
(bonded  and  torsional  terms  plus  internal  vdW)  produces  a 
small  contribution  that  favors  the  ‘out’  macrostate.  The 
dominant  terms  are  the  electrostatic  compensation  of 
A Wgie  vs.  AGS  ele,  and  conformational  entropy.  The  latter 
reveals  a  topology  of  a  wider  ‘out’  basin  and  allows  for 
greater  dynamic  connectivity  among  microstates. 

The  agreement  obtained  from  the  LRA  model  with  the 
experimental  estimate  is  encouraging  from  the  perspec¬ 
tive  of  not  having  to  scale  a  in  eq.  (7)  from  the  initial  value 
of  0.5.  Generally,  for  the  application  of  LRA,  the  a  priori 
determination  of  a  clearly  hinders  its  use,  although  if  data 
are  available,  mutational  free  energies  or  ligand-binding 
affinities  can  be  used  to  fit  both  scaling  parameters  of  eq. 
(7). 

The  next  computational  approach  is  the  rescoring  of  the 
MD  trajectories  by  removing  the  explicit  water  and  calcu¬ 
lating  the  solvent  polarization  term  with  either  the  GB 
model  or  the  Poisson  equation.  Table  I  shows  that  values  of 
AGS  eie  for  both  implicit  solvent  models  are  underesti¬ 
mated  in  comparison  with  the  LRA  calculation.  This  is 
most  evident  for  the  more  rigorous  Poisson  equation.  The 
observed  reduction  in  differentiation  of  solvent  polariza¬ 
tion  between  the  macrosates  is  the  outcome  of  transition¬ 
ing  conformations  that  were  thermally  equilibrated  on  an 
explicit  solvent  potential  energy  surface  to  a  non-optimal 
implicit  solvent  environment.  The  solvent  perturbation 
incorrectly  alters  the  energy  landscape  without  any  confor¬ 
mational  relaxation  to  adapt  to  the  new  topology.  In  other 
words,  conformations  culled  at  298K  in  explicit  solvent  are 
positioned  at  a  higher  energy  on  the  potential  energy 


surface  from  changing  to  implicit  solvent.  Errors  arise 
principally  from  charged  residues  (e.g.,  Arg-289),  and  the 
lack  of  cancellation  of  relative  error  in  AGS  ele  is  due  to 
conformational  dependence  of  the  electric  field.  Reconcilia¬ 
tion  of  the  implicit  solvent  results  by  simple  scaling  the 
homogeneous  ep  to  account  for  dipolar  reorientation  will 
fail  to  provide  consistent  improvement.  As  commented 
below,  the  calculations  illustrate  the  difficulty  of  obtaining 
accurate  and  reliable  results  from  the  approach  of  rescor¬ 
ing  conformations  with  a  solvent  model  different  from  the 
one  used  to  generate  the  ensemble. 

The  third  set  of  calculations  is  the  MD-GB/SA  model 
simulations,  which  show  accuracy  in  A Gmacro  comparable 
with  that  of  the  LRA.  This  result  is  promising,  yet  inspec¬ 
tion  of  the  electrostatic  free  energy  terms  reveals  signifi¬ 
cant  differences  in  sampling  microstates.  Conformations 
found  from  exploring  the  GB  potential  surface  are  opti¬ 
mized  as  a  result  of  instantaneous  solvent  dipolar  relax¬ 
ation  embedded  in  the  mean-field  approximation,  while 
explicit  water  calculations  treat  dipolar  reorientation. 
This  distinction  is  apparent  from  scatter  plots  of  FDP  vs. 
GB  (Fig.  2).  Scoring  conformations  by  FDP  reveals  roughly 
the  same  range  of  values  as  the  ensemble  generated  from 
MD-GB/SA  and  from  explicit  water  calculations,  while 
scoring  the  MD-GB/SA  ensemble  by  GB  favors  a  shift  on 
the  order  of  100  kcal/mol  from  the  ideal  correlation  be¬ 
tween  GB  and  FDP.  The  latter  is  inconsistent  with  the 
explicit  water  results;  the  MD-GB/SA  model  simply 
searches  for  conformations  optimized  by  an  implicit  reac¬ 
tion  field.  The  strength  of  interaction  between  polarization 
charges  of  the  protein  and  solvent  scales  as  (l/ep  -  1/80) 
with  the  maximum  value  obtained  for  e  =  1  and  charged 
residues  in  their  ionized  state.  The  indirect  success  of  GB 
simulations  without  scaling  the  protein  dielectric  constant 
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Fig.  2.  Scatter  plots  of  FDP  vs.  GB  solvation  free  energies  calculated  from  scoring  molecular  dynamics  simulation  trajectories  using  different  solvent 
models.  The  basin  of  ‘in’  conformations  is  shown  in  blue  and  the  ‘out’  in  red.  (a)  Explicit  solvent  results;  correlation  coefficient  (t2)  is  0.9  for  the  ‘in’  basin 
and  0.6  for  the  ‘out.’  (b)  MD-GB/SA  model;  i2  =  0.8  for  both  basins,  (c)  ASP  solvent  model;  r2  =  0.9  and  0.8  for  the  ‘in’  and  ‘out’  conformations, 
respectively,  (d)  The  s  =  r  solvent  model;  r2  =  0.8  and  0.7  for  the  ‘in’  and  ‘out’  conformations,  respectively.  Each  plot  shows  an  ideal  correlation  and  linear 
fits  to  the  computed  data. 


is  due  to  sufficient  configurational  averaging  of  the  basins 
and  their  smooth-funnel  topology.18 

While  the  manifold  of  microstates  sampled  by  the  ex¬ 
plicit  solvent  model  and  MD-GB/SA  are  different,  similar 
results  in  A Gmacro  are  predominantly  due  to  the  effect  of 
the  net  electrostatic  component  on  conformational  sam¬ 
pling.  From  the  explicit  water  simulation,  A  Wele  4-  AGS  ele  = 
14.6  kcal/mol,  while  from  the  MD-GB/SA,  the  sum  equals 
11.5  kcal/mol.  The  two  simulation  models  show  similar 
trends,  whereas  models  of  rescoring  conformations  either 
underestimate  or  failed  to  properly  offset  Coulombic  charg¬ 
ing  by  solvation.  Again  this  supports  the  notion  of  obtain¬ 
ing  self-consistency  between  generating  ensembles  and 
scoring  them. 

The  overall  scatter  from  a  perfect  linear  fit  of  FDP  vs. 
GB  reflects  insufficient  parameterization  of  earlier  GB 
models23  for  protein  simulations  and  can  be  improved  with 
more  recent  models.28  For  the  explicit  solvent  calculations, 
the  correlation  coefficient  is  0.9  for  the  ‘in’  macrostate  and 
0.6  for  the  ‘out’,  while  GB  simulation  yields  0.8  for  both 
macrostates. 

From  Table  I,  the  ASP  model  and  the  distance-depen¬ 
dent  screening  of  charges  are  unsuccessful  at  producing 


accurate  results.  Both  models  show  significant  protein 
distortions  reflected  in  AEnp  and  are  no  more  reliable  than 
vacuum  calculations.  Other  modeling  studies  have  noted 
protein  distortions  that  result  from  using  similar  solvent 
models.29  Despite  the  excellent  agreement  between  GB 
and  FDP  for  revaluating  the  trajectories  generated  by  the 
ASP  model  (Fig.  2),  the  calculations  show  that  the  ‘out’ 
basin  is  more  favorable  by  >60  kcal/mol.  In  comparison, 
dielectric  screening  correctly  discriminates  the  ‘in’  basin 
by  GB  evaluation,  yet  consistent  scoring  by  the  function 
used  to  generate  the  ensemble  fails.  Incorrect  modeling  of 
solvent  polarization  and  conformational  distortions  ren¬ 
ders  both  models  poor  choices  for  sampling  microstates. 
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